knitr::opts_chunk$set(echo = T)
library(catalogueR)
library(dplyr)
root <- getwd()
hpc_root <- "/sc/arion/projects/pd-omics/brian/COVID19_GWAS"DAT <- readxl::read_excel("data/GCST90000255and6_GRCh38_Plessthan1eMinus5_combined_pruned.xlsx") %>%
dplyr::rename(SNP=`SNP ID`,
CHR=chromosome,
POS=base_pair_location,
A1=effect_allele,
A2=other_allele,
Effect=beta,
StdErr=standard_error,
P=p_value) %>%
dplyr::mutate(Locus=SNP)
createDT(DAT)Since catalogueR expects a list of files with per-Locus summary stats, we’ll just save each SNP as its own “locus” file.
Annotate tissues using metadata
Count the number of unique SNPs, eQTL genes, significant eGenes per tissue/cell-type.
Note: eQTL Catalogue only provides raw p-values, so you’ll need to infer your own signficance correction threshold. For simplicity, I use p < 5e-8 here.
GWAS.QTL <- dplyr::mutate(GWAS.QTL, FDR.QTL=stats::p.adjust(p = pvalue.QTL, method = "fdr"))
# p_thresh <- 5e-8
fdr_thresh <- 0.05
counts <- GWAS.QTL %>%
dplyr::group_by(System, Tissue_group, Tissue) %>%
dplyr::summarise(variants=dplyr::n_distinct(SNP),
eVariants=dplyr::n_distinct(SNP[FDR.QTL<fdr_thresh]),
genes=dplyr::n_distinct(gene.QTL),
eGenes=dplyr::n_distinct(gene.QTL[FDR.QTL<fdr_thresh]),
eGenes_names=paste(unique(gene.QTL[FDR.QTL<fdr_thresh]), collapse = ", "),
studies=dplyr::n_distinct(Study),
studies_names=paste(unique(Study), collapse = ", ")) %>%
dplyr::arrange(desc(variants),desc(eVariants),
desc(genes), desc(eGenes))## `summarise()` regrouping output by 'System', 'Tissue_group' (override with `.groups` argument)
library(ggplot2)
GWAS.QTL$Tissue <- factor(GWAS.QTL$Tissue, levels = unique(counts$Tissue), ordered = T)
gp <- ggplot(data = subset(GWAS.QTL, FDR.QTL< fdr_thresh),
aes(x= gene.QTL, y= Tissue, color=P, size= FDR.QTL)) +
geom_point(alpha=1) +
scale_color_viridis_c(direction = -1) +
facet_grid(facets = System ~.,
space = "free_y",
scales = "free_y") +
guides(size = guide_legend(reverse=T)) +
labs(title = "Variant overlap: COVID-19 GWAS vs. eQTL Catalogue",
subtitle = paste("QTL FDR <",fdr_thresh)) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill="grey20"),
strip.text = element_text(color="white"))
print(gp)## Saving 9 x 7.5 in image